# Install and load the readxl package if not already installed
if (!requireNamespace("readxl", quietly = TRUE)) {
install.packages("readxl")
}
library(readxl)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# Read the Excel file and convert it to a dataframe
df <- read_excel("Hornets data.xlsx")
## Warning: Expecting logical in AN6909 / R6909C40: got
## 'CA_AB_Webb_1959_08_02_001'
# Filter the dataframe and remove rows with missing longitude or latitude
target_df <- df %>%
filter(stateProvince %in% c('British Columbia', 'Canada - British Columbia (BC)')) %>%
filter(!is.na(decimalLongitude) | !is.na(decimalLatitude))
if (!requireNamespace("maptools", quietly = TRUE)) {
install.packages("maptools", repos = "https://packagemanager.posit.co/cran/2023-10-13")
}
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
## Checking rgeos availability: FALSE
# Load required packages
library(sp)
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 3.2-9
## Loading required package: spatstat.random
## spatstat.random 3.2-3
## Loading required package: spatstat.explore
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## spatstat.explore 3.2-7
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.2-11
## Loading required package: spatstat.linnet
## spatstat.linnet 3.1-5
##
## spatstat 3.0-8
## For an introduction to spatstat, type 'beginner'
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(maptools)
##
## Attaching package: 'maptools'
## The following object is masked from 'package:sp':
##
## sp2Mondrian
# Import BC_Covariates
load("BC_Covariates.Rda")
window = DATA$Window
library(ggplot2)
# Plot longitude and latitude coordinates
ggplot(target_df, aes(x = decimalLongitude, y = decimalLatitude)) +
geom_point(color = "blue") + # Add points
labs(x = "Longitude", y = "Latitude", # Label axes
title = "Location Coordinates for Hornet within BC", # Add title
color = "Data Points") + # Legend title
scale_color_manual(values = "blue", name = "Data Points") + # Legend color and title
theme_minimal() + # Apply minimal theme
theme(
panel.grid = element_blank(), # Remove grid lines
axis.line = element_line(color = "black") # Color axis lines
)
library(sp)
# Define the projection string
proj_string <- "+proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs"
# Create a SpatialPoints object from latitude and longitude
points <- SpatialPoints(coords = target_df[, c("decimalLongitude", "decimalLatitude")],
proj4string = CRS("+proj=longlat +datum=WGS84"))
# Transform the points to BC Albers projection
points_bc_albers <- spTransform(points, CRS(proj_string))
# Extract X and Y coordinates
X <- coordinates(points_bc_albers)[, 1]
Y <- coordinates(points_bc_albers)[, 2]
# Create a dataframe with X and Y columns
df_bc_albers <- data.frame(X = X, Y = Y)
# Now, 'df_bc_albers' contains the X and Y columns in the BC Albers projection
# remove outside points
valid_points <- df_bc_albers[inside.owin(df_bc_albers$X, df_bc_albers$Y, as.owin(DATA$Window)), ]
hornet_ppp <- ppp(
x = valid_points$X, # X coordinates of valid points
y = valid_points$Y, # Y coordinates of valid points
window = as.owin(DATA$Window) # Observation window
)
## Warning: data contain duplicated points
plot(hornet_ppp)
From the plot, we can see most hornets tend to gather in lower BC, showing high intensity especially at the southwest corner. But for other part, especially for upper BC, there is no evidence of correlation in hornet location from eyes, indicating the distribution is inhomogeneous.
win_km <- rescale(Window(hornet_ppp), 1000, "km")
# Intensity in trees/km^2
npoints(hornet_ppp)/area(win_km)
## [1] 0.001371983
I dont think the estimated intensity trustworthy. Because from the plot we cannot tell the hornet location is homogeneous.
Q <- quadratcount(hornet_ppp,
nx = 10,
ny = 10)
quadrat.test(Q)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 28479, df = 63, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 64 tiles (irregular windows)
The small p-value suggests that there is a significant deviation from homogeneity.
plot(Q, main = "Quadrat Plot with Hornet Points", xlab = "X Coordinate", ylab = "Y Coordinate")
points(hornet_ppp, pch = 20, col = "green") # Overlay hornet points
intensity(Q)
## tile
## Tile row 1, col 1 Tile row 1, col 2 Tile row 1, col 3 Tile row 1, col 4
## 0.000000e+00 1.734704e-10 0.000000e+00 2.561876e-10
## Tile row 1, col 5 Tile row 1, col 6 Tile row 1, col 7 Tile row 2, col 2
## 6.279768e-10 0.000000e+00 0.000000e+00 0.000000e+00
## Tile row 2, col 3 Tile row 2, col 4 Tile row 2, col 5 Tile row 2, col 6
## 0.000000e+00 0.000000e+00 0.000000e+00 1.374831e-10
## Tile row 2, col 7 Tile row 3, col 2 Tile row 3, col 3 Tile row 3, col 4
## 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## Tile row 3, col 5 Tile row 3, col 6 Tile row 3, col 7 Tile row 4, col 3
## 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## Tile row 4, col 4 Tile row 4, col 5 Tile row 4, col 6 Tile row 4, col 7
## 4.582771e-11 0.000000e+00 0.000000e+00 3.146336e-10
## Tile row 5, col 3 Tile row 5, col 4 Tile row 5, col 5 Tile row 5, col 6
## 0.000000e+00 2.774770e-10 5.499325e-10 4.582771e-11
## Tile row 5, col 7 Tile row 6, col 2 Tile row 6, col 3 Tile row 6, col 4
## 0.000000e+00 0.000000e+00 0.000000e+00 4.878190e-11
## Tile row 6, col 5 Tile row 6, col 6 Tile row 6, col 7 Tile row 6, col 8
## 9.165542e-11 1.237348e-09 4.591895e-11 0.000000e+00
## Tile row 7, col 2 Tile row 7, col 3 Tile row 7, col 4 Tile row 7, col 5
## 0.000000e+00 0.000000e+00 7.697947e-11 4.666335e-11
## Tile row 7, col 6 Tile row 7, col 7 Tile row 7, col 8 Tile row 7, col 9
## 9.165542e-11 9.165542e-11 5.162269e-11 0.000000e+00
## Tile row 8, col 4 Tile row 8, col 5 Tile row 8, col 6 Tile row 8, col 7
## 0.000000e+00 0.000000e+00 9.195917e-11 7.332433e-10
## Tile row 8, col 8 Tile row 8, col 9 Tile row 8, col 10 Tile row 9, col 4
## 3.712044e-09 1.500022e-09 2.603879e-10 0.000000e+00
## Tile row 9, col 5 Tile row 9, col 6 Tile row 9, col 7 Tile row 9, col 8
## 1.741078e-09 2.327311e-09 5.499325e-10 4.491115e-09
## Tile row 9, col 9 Tile row 9, col 10 Tile row 10, col 5 Tile row 10, col 6
## 1.420659e-09 7.075868e-10 1.130684e-09 5.472572e-08
## Tile row 10, col 7 Tile row 10, col 8 Tile row 10, col 9 Tile row 10, col 10
## 1.855185e-08 2.875371e-09 2.577527e-09 0.000000e+00
plot ( intensity ( Q , image = T))
lambda_u <- mean (Q)
R_Poisson <- rpois (n = length (Q) ,
lambda = lambda_u)
hist(Q)
hist(R_Poisson)
These data are clearly inhomogeneous.
plot(density(hornet_ppp, sigma = bw.ppl), # Likelihood Cross Validation Bandwidth Selection
ribbon = F,
main = "Hornet Kernel Estimation")
R <- bw.ppl(hornet_ppp)
#Calculate test statistic
LR <- scanLRTS(hornet_ppp, r = R)
#Plot the output
plot(LR,main='Hotspot Analysis')
plot(window,add = TRUE)
#Estimate a strictly positive density
lambda_hornet_pos <- density(hornet_ppp,
sigma=bw.ppl,
positive=TRUE)
#Simulation envelope (with points drawn from the estimated intensity)
E_hornet_inhom <- envelope(hornet_ppp,
Kinhom,
simulate = expression(rpoispp(lambda_hornet_pos)),
correction="border",
rank = 1,
nsim = 19,
fix.n = TRUE)
## Warning in envelope.ppp(hornet_ppp, Kinhom, simulate =
## expression(rpoispp(lambda_hornet_pos)), : fix.n and fix.marks were ignored,
## because 'simulate' was given
## Generating 19 simulations by evaluating expression ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
# visualise the results
par(mfrow = c(1,2))
plot(E_hornet_inhom,
main = "",
lwd = 2)
# Zoom in on range where significant deviations appear
plot(E_hornet_inhom,
xlim = c(0,30000),
main = "",
lwd = 2)
When corrected for inhomogeneity, significant clustering only appears to exist in and around 0-22000 meters. For longer distance, we cannot indicate a significant evidence of correlation in hornet distribution.
pcf_hornet_inhom <- envelope(hornet_ppp,
pcfinhom,
simulate = expression(rpoispp(lambda_hornet_pos)),
rank = 1,
nsim = 19)
## Generating 19 simulations by evaluating expression ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
par(mfrow = c(1,3))
plot(pcf_hornet_inhom)
# Zoom in on range where significant deviations appear
plot(pcf_hornet_inhom,
xlim = c(0,10000),
main = "",
lwd = 2)
plot(pcf_hornet_inhom,
xlim = c(15000,50000),
main = "",
lwd = 2)
There appear to be more hornet than expected by random chance between ∼ 0 - 6000 meters and less hornet between 15000 - 50000 meters. Except that, the locations of hornet appear not to exhibit any significant correlations.
plot(DATA$Elevation, main = "Elevation Map of British Columbia", xlab = "Longitude", ylab = "Latitude")
# Plot park locations on top of elevation map
plot(hornet_ppp, add = TRUE,cex = 0.4)
# Define elevation classes based on quantiles
elevation_classes <- cut(DATA$Elevation, breaks = 5, labels = FALSE)
# Plot the elevation class image
image(elevation_classes, col = terrain.colors(5), main = "Elevation Classes")
plot(hornet_ppp, add = TRUE,cex = 0.4)
From plot, we could find that most hornets gather in Class 1 and 2.
# Median elevation within British Columbia
median_bc_elevation <- median(DATA$Elevation, na.rm = TRUE)
# Median elevation at park locations
median_hornet_elevation <- median(DATA$Elevation[hornet_ppp], na.rm = TRUE)
# Print results
print(paste("Median elevation within British Columbia:", median_bc_elevation))
## [1] "Median elevation within British Columbia: 1096.85723142938"
print(paste("Median elevation at hornet locations:", median_hornet_elevation))
## [1] "Median elevation at hornet locations: 74.4980033424321"
BC_elevation_density <- density(as.data.frame(DATA$Elevation, na.rm = TRUE)$value)
# Extract elevation values from DATA$Elevation using hornet_ppp
elevation_values <- as.vector(DATA$Elevation[hornet_ppp])
# Generate a kernel density estimate of the distribution of elevation values within hornet
hornet_elevation_density <- density(elevation_values, na.rm = TRUE)
plot(BC_elevation_density, col = "blue", main = "Kernel Density Estimate of Elevation in British Columbia", xlab = "Elevation",ylim = c(0, max(hornet_elevation_density$y)))
# Overlay the kernel density estimate of elevation values for park locations
lines(hornet_elevation_density, col = "red")
# Add legend
legend("topright", legend = c("British Columbia", "hornet Locations"), col = c("blue", "red"), lty = 1, cex = 0.8)
plot(DATA$Forest, main = "Forest Map of British Columbia", xlab = "Longitude", ylab = "Latitude")
# Plot park locations on top of elevation map
plot(hornet_ppp, add = TRUE,cex = 0.4)
# Define elevation classes based on quantiles
forest_classes <- cut(DATA$Forest, breaks = 5, labels = FALSE)
# Plot the elevation class image
image(forest_classes, col = terrain.colors(5), main = "Forest Classes")
plot(hornet_ppp, add = TRUE,cex = 0.4)
From plot, we could find that most hornets gather in Class 4 and 5.
# Median forestdensity within British Columbia
median_bc_forest <- median(DATA$Forest, na.rm = TRUE)
# Median forest density at park locations
median_hornet_forest <- median(DATA$Forest[hornet_ppp], na.rm = TRUE)
# Print results
print(paste("Median Forest Density within British Columbia:", median_bc_forest))
## [1] "Median Forest Density within British Columbia: 50.1524677276611"
print(paste("Median Forest Density at hornet locations:", median_hornet_forest))
## [1] "Median Forest Density at hornet locations: 19.1723461151123"
BC_forest_density <- density(as.data.frame(DATA$Forest, na.rm = TRUE)$value)
# Extract forest values from DATA$Forest using hornet_ppp
forest_values <- as.vector(DATA$Forest[hornet_ppp])
# Generate a kernel density estimate of the distribution of forest values within hornet
hornet_forest_density <- density(forest_values, na.rm = TRUE)
plot(BC_forest_density, col = "blue", main = "Kernel Density Estimate of Forest in British Columbia", xlab = "Forest",ylim = c(0, max(BC_forest_density$y)))
# Overlay the kernel density estimate of forest values for hornet locations
lines(hornet_forest_density, col = "red")
# Add legend
legend("topright", legend = c("British Columbia", "hornet Locations"), col = c("blue", "red"), lty = 1, cex = 0.8)
### Distance to Water Map
plot(DATA$Dist_Water, main = "Distance to Water Map of British Columbia", xlab = "Longitude", ylab = "Latitude")
# Plot park locations on top of elevation map
plot(hornet_ppp, add = TRUE,cex = 0.4,col = 'white')
From the graph we can find that the distance is almost pretty low from
water within BC, so i tend not to go deeper into the single
indicator.
plot(DATA$HFI, main = "HFI Map of British Columbia", xlab = "Longitude", ylab = "Latitude")
# Plot park locations on top of elevation map
plot(hornet_ppp, add = TRUE,cex = 0.4, col='white')
From the graph we can find that HFI is almost pretty low within BC, so i
tend not to go deeper into the single indicator.
elev <- DATA$Elevation
hornet_elev_rho <- rhohat(hornet_ppp, elev)
plot(hornet_elev_rho,xlim = c(0,max(elev)))
fors <- DATA$Forest
hornet_forest_rho <- rhohat(hornet_ppp, fors)
plot(hornet_forest_rho,xlim = c(0,max(fors)))
dist_water <- DATA$Dist_Water
hornet_dist_rho <- rhohat(hornet_ppp, dist_water)
plot(hornet_dist_rho,xlim = c(0,max(dist_water)))
hfi <- DATA$HFI
So we could find that the three indicators cannot be linear with hornet distribution, which lead to the formula below.
fit <- ppm(hornet_ppp ~ elev + fors + dist_water + hfi , data = DATA)
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of
## 4230) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## hornet_ppp, trend = ~elev + fors + dist_water + hfi,
fit
## Nonstationary Poisson process
## Fitted to point pattern dataset 'hornet_ppp'
##
## Log intensity: ~elev + fors + dist_water + hfi
##
## Fitted trend coefficients:
## (Intercept) elev fors dist_water hfi
## -2.085013e+01 -2.136346e-03 4.101943e-04 -6.229289e-05 6.388678e+00
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.085013e+01 1.394446e-01 -2.112344e+01 -2.057683e+01 ***
## elev -2.136346e-03 1.018394e-04 -2.335948e-03 -1.936745e-03 ***
## fors 4.101943e-04 1.347051e-03 -2.229978e-03 3.050367e-03
## dist_water -6.229289e-05 1.424637e-05 -9.021527e-05 -3.437051e-05 ***
## hfi 6.388678e+00 1.516429e-01 6.091463e+00 6.685893e+00 ***
## Zval
## (Intercept) -149.5226796
## elev -20.9775976
## fors 0.3045127
## dist_water -4.3725438
## hfi 42.1297487
## Problem:
## Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of 4230)
## of the quadrature points
From the Ztest, we’d better remove forest
fit_1 <- ppm(hornet_ppp ~ elev + dist_water + hfi, data = DATA)
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of
## 4230) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## hornet_ppp, trend = ~elev + dist_water + hfi, data = list(
fit_1
## Nonstationary Poisson process
## Fitted to point pattern dataset 'hornet_ppp'
##
## Log intensity: ~elev + dist_water + hfi
##
## Fitted trend coefficients:
## (Intercept) elev dist_water hfi
## -2.082381e+01 -2.135456e-03 -6.280277e-05 6.366200e+00
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.082381e+01 1.092606e-01 -2.103796e+01 -2.060967e+01 ***
## elev -2.135456e-03 1.018300e-04 -2.335039e-03 -1.935872e-03 ***
## dist_water -6.280277e-05 1.414419e-05 -9.052487e-05 -3.508066e-05 ***
## hfi 6.366200e+00 1.323565e-01 6.106786e+00 6.625614e+00 ***
## Zval
## (Intercept) -190.588484
## elev -20.970799
## dist_water -4.440181
## hfi 48.098883
## Problem:
## Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of 4230)
## of the quadrature points
fit_2 <- ppm(hornet_ppp ~ elev + I(elev^2)+ dist_water + hfi + I(hfi^2), data = DATA)
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of
## 4230) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## hornet_ppp, trend = ~elev + I(elev^2) + dist_water +
fit_2
## Nonstationary Poisson process
## Fitted to point pattern dataset 'hornet_ppp'
##
## Log intensity: ~elev + I(elev^2) + dist_water + hfi + I(hfi^2)
##
## Fitted trend coefficients:
## (Intercept) elev I(elev^2) dist_water hfi
## -2.120567e+01 -3.327566e-03 9.273141e-07 -5.705994e-05 1.034704e+01
## I(hfi^2)
## -4.270298e+00
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.120567e+01 1.406367e-01 -2.148132e+01 -2.093003e+01 ***
## elev -3.327566e-03 2.021561e-04 -3.723784e-03 -2.931347e-03 ***
## I(elev^2) 9.273141e-07 1.358787e-07 6.609967e-07 1.193631e-06 ***
## dist_water -5.705994e-05 1.387989e-05 -8.426403e-05 -2.985586e-05 ***
## hfi 1.034704e+01 5.855823e-01 9.199316e+00 1.149476e+01 ***
## I(hfi^2) -4.270298e+00 5.922054e-01 -5.430999e+00 -3.109596e+00 ***
## Zval
## (Intercept) -150.783353
## elev -16.460374
## I(elev^2) 6.824573
## dist_water -4.110979
## hfi 17.669652
## I(hfi^2) -7.210838
## Problem:
## Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of 4230)
## of the quadrature points
The coefficients are all statistically significant.
#Plot the model predictions
plot(fit_2,
se = FALSE,
superimpose = FALSE)
#Overlay the B. pendula locations
plot(hornet_ppp,
pch = 16,
cex = 0.4,
cols = "white",
add = TRUE)
plot(hornet_ppp,
pch = 16,
cex = 0.3,
cols = "black",
add = TRUE)
AIC(fit_2);AIC(fit_1)
## [1] 47779.17
## [1] 47850.98
BIC(fit_2); BIC(fit_1)
## [1] 47810.2
## [1] 47871.66
We could find that AIC and BIC fall down when we remove one feature. So the extra complexity is not well supported by the data.
#Run the quadrat test
quadrat.test(fit_2, nx = 10, ny = 10)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of fitted Poisson model 'fit_2' using quadrat counts
##
## data: data from fit_2
## X2 = 557.63, df = 58, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 64 tiles (irregular windows)
The small p value tells us that there’s a significant deviation from our model’s predictions. While this is useful for suggesting that our model has room for improvement, it provides us with no direction on how to do so (e.g., missing parameters, model mispecification (e.g., polynomial vs. linear), a lack of independence, non-stationarity, etc…).
#Calculate the residuals
res <- residuals(fit_2)
## Warning: Some infinite, NA or NaN increments were removed
na_indexes <- which(is.na(res$val))
# If you need to exclude NA values explicitly
if (length(na_indexes) > 0) {
res <- res[-na_indexes]
plot(res, cols='transparent')
} else {
plot(res, cols='transparent')
}
We can see some pattern in the residual plot. We still need to check
#Calculate the partial residuals as a function of elevation
par_res_elev0 <- parres(fit_2, "elev")
## Warning: Some infinite, NA or NaN increments were removed
#Calculate the relative intensity as a function of for
par_res_hfi0 <- parres(fit_2, "hfi")
## Warning: Some infinite, NA or NaN increments were removed
## Warning: Values for 2 query points lying outside the pixel image domain were
## estimated by projection to the nearest pixel
par_res_water0 <- parres(fit_2,"dist_water")
## Warning: Some infinite, NA or NaN increments were removed
#Side by side plotting
par(mfrow = c(1,3))
plot(par_res_elev0,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Elevation (m)")
plot(par_res_hfi0,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Forest")
plot(par_res_water0,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Distance to Water")
The plot shows that each feature could explain the intensity well.
library(splines)
#Fit the PPP model
fit_smooth <- ppm(hornet_ppp ~ bs(elev,4) + bs(hfi, 7) + bs(dist_water,2), data = DATA, use.gam = TRUE)
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of
## 4230) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## hornet_ppp, trend = ~bs(elev, 4) + bs(hfi, 7) + bs(dist_water,
## Warning in bs(dist_water, 2): 'df' was too small; have used 3
## Warning in bs(dist_water, 2): 'df' was too small; have used 3
fit_smooth
## Nonstationary Poisson process
## Fitted to point pattern dataset 'hornet_ppp'
##
## Log intensity: ~bs(elev, 4) + bs(hfi, 7) + bs(dist_water, 2)
##
## Fitted trend coefficients:
## (Intercept) bs(elev, 4)1 bs(elev, 4)2 bs(elev, 4)3
## -22.6986660 -2.0545674 -2.4037646 -4.0002691
## bs(elev, 4)4 bs(hfi, 7)1 bs(hfi, 7)2 bs(hfi, 7)3
## -11.3093207 0.9351316 2.3075726 2.4558059
## bs(hfi, 7)4 bs(hfi, 7)5 bs(hfi, 7)6 bs(hfi, 7)7
## 5.4655838 6.7798595 6.9551590 9.2755929
## bs(dist_water, 2)1 bs(dist_water, 2)2 bs(dist_water, 2)3
## 0.7654902 -3.0664677 -2.3735535
##
## For standard errors, type coef(summary(x))
## Problem:
## Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of 4230)
## of the quadrature points
coef(summary(fit_smooth))
## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -22.6986660 0.7302310 -24.12989251 -21.2674395 ***
## bs(elev, 4)1 -2.0545674 0.2503921 -2.54532701 -1.5638079 ***
## bs(elev, 4)2 -2.4037646 0.5240757 -3.43093409 -1.3765952 ***
## bs(elev, 4)3 -4.0002691 1.2726075 -6.49453398 -1.5060043 **
## bs(elev, 4)4 -11.3093207 5.6102003 -22.30511111 -0.3135302 *
## bs(hfi, 7)1 0.9351316 0.9798759 -0.98538988 2.8556531
## bs(hfi, 7)2 2.3075726 0.8685392 0.60526719 4.0098781 **
## bs(hfi, 7)3 2.4558059 0.7503720 0.98510380 3.9265081 **
## bs(hfi, 7)4 5.4655838 0.7615734 3.97292739 6.9582403 ***
## bs(hfi, 7)5 6.7798595 0.7676102 5.27537107 8.2843479 ***
## bs(hfi, 7)6 6.9551590 0.7374069 5.50986799 8.4004500 ***
## bs(hfi, 7)7 9.2755929 0.7572240 7.79146117 10.7597246 ***
## bs(dist_water, 2)1 0.7654902 0.3477852 0.08384376 1.4471366 *
## bs(dist_water, 2)2 -3.0664677 1.0435042 -5.11169839 -1.0212370 **
## bs(dist_water, 2)3 -2.3735535 2.0392890 -6.37048648 1.6233794
## Zval
## (Intercept) -31.0842253
## bs(elev, 4)1 -8.2053990
## bs(elev, 4)2 -4.5866746
## bs(elev, 4)3 -3.1433644
## bs(elev, 4)4 -2.0158497
## bs(hfi, 7)1 0.9543368
## bs(hfi, 7)2 2.6568435
## bs(hfi, 7)3 3.2727845
## bs(hfi, 7)4 7.1767000
## bs(hfi, 7)5 8.8324246
## bs(hfi, 7)6 9.4319143
## bs(hfi, 7)7 12.2494706
## bs(dist_water, 2)1 2.2010432
## bs(dist_water, 2)2 -2.9386251
## bs(dist_water, 2)3 -1.1639123
#Calculate the partial residuals as a function of elevation
par_res_elev <- parres(fit_smooth, "elev")
## Warning: Some infinite, NA or NaN increments were removed
## Warning in bs(hfi, 7): all interior knots match left boundary knot
## Warning in bs(hfi, 7): all interior knots match right boundary knot
## Warning in bs(dist_water, 2): 'df' was too small; have used 3
#Calculate the relative intensity as a function of for
par_res_hfi <- parres(fit_smooth, "hfi")
## Warning: Some infinite, NA or NaN increments were removed
## Warning: Values for 2 query points lying outside the pixel image domain were
## estimated by projection to the nearest pixel
## Warning in bs(elev, 4): all interior knots match left boundary knot
## Warning in bs(elev, 4): all interior knots match right boundary knot
## Warning in bs(dist_water, 2): 'df' was too small; have used 3
par_res_water <- parres(fit_smooth,"dist_water")
## Warning: Some infinite, NA or NaN increments were removed
## Warning in bs(elev, 4): all interior knots match left boundary knot
## Warning in bs(elev, 4): all interior knots match right boundary knot
## Warning in bs(hfi, 7): all interior knots match left boundary knot
## Warning in bs(hfi, 7): all interior knots match right boundary knot
## Warning in bs(dist_water, 2): 'df' was too small; have used 3
#Side by side plotting
par(mfrow = c(1,3))
plot(par_res_elev,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Elevation (m)")
plot(par_res_hfi,
legend = FALSE,
lwd = 2,
main = "",
xlab = "HFI")
plot(par_res_water,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Distance to Water")
It is not good, but better than previous one.
AIC(fit_2);AIC(fit_smooth)
## [1] 47779.17
## [1] 47645
BIC(fit_2);BIC(fit_smooth)
## [1] 47810.2
## [1] 47722.56
anova(fit_2, fit_smooth, test = "LRT")
## Warning: Models were re-fitted with use.gam=TRUE
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of
## 4230) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## hornet_ppp, trend = ~elev + I(elev^2) + dist_water +
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.31% (13 out of
## 4230) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## hornet_ppp, trend = ~bs(elev, 4) + bs(hfi, 7) + bs(dist_water,
## Warning in bs(dist_water, 2): 'df' was too small; have used 3
## Warning in bs(dist_water, 2): 'df' was too small; have used 3
## Analysis of Deviance Table
##
## Model 1: ~elev + I(elev^2) + dist_water + hfi + I(hfi^2) Poisson
## Model 2: ~bs(elev, 4) + bs(hfi, 7) + bs(dist_water, 2) Poisson
## Npar Df Deviance Pr(>Chi)
## 1 6
## 2 15 9 152.18 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The ANOVA results indicate that Model 2 has significantly better fit to the data than Model 1, as evidenced by the large difference in deviance (516.81) and the highly significant p-value (< 2.2e-16), which is well below any conventional significance level. All lines of evidence point towards these more complex models being a better fit to the data. The final model is too long to write down, but we can visualise the predictions just as before.
#Calculate the residuals
res1 <- residuals(fit_smooth)
## Warning: Some infinite, NA or NaN increments were removed
na_indexes1 <- which(is.na(res1$val))
# If you need to exclude NA values explicitly
if (length(na_indexes1) > 0) {
res1 <- res[-na_indexes1]
plot(res1, cols='transparent')
} else {
plot(res1, cols='transparent')
}
#Plot the model predictions
plot(fit_smooth,
se = FALSE,
superimpose = FALSE,
main = "Estimated Hornet intensity")
plot(hornet_ppp,
pch = 16,
cex = 0.6,
cols = "white",
add = TRUE)
plot(hornet_ppp,
pch = 16,
cex = 0.5,
cols = "black",
add = TRUE)